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Multi-reference calculations along the lines of the Generator Coordinate Method or the restoration 
of broken symmetries within the nuclear Energy Density Functional (EDF) framework are becoming 
a standard tool in nuclear structure physics. These calculations rely on the extension of a single- 
reference energy functional, of the Gogny or the Skyrme types, to non-diagonal energy kernels. 
There is no rigorous constructive framework for this extension so far. The commonly accepted 
way proceeds by formal analogy with the expressions obtained when applying the generalized Wick 
theorem to the non-diagonal matrix element of a Hamilton operator between two product states. It is 
pointed out that this procedure is ill-defined when extended to EDF calculations as the generalized 
Wick theorem is taken outside of its range of applicability. In particular, such a procedure is 
responsible for the appearance of spurious divergences and steps in multi-reference EDF energies, 
as was recently observed in calculations restoring particle number or angular momentum. In the 
present work, we give a formal analysis of the origin of this problem for calculations with and without 
pairing, i.e. constructing the density matrices from either Slater determinants or quasi-particle 
vacua. We propose a method to regularize non-diagonal energy kernels such that divergences and 
steps are removed from multi-reference EDF energies. Such a removal is a priori quasi-particle-basis 
dependent. A special feature of the method we use to proceed to the actual regularization is that it 
singles out one basis among all possible ones. The regularization method is applicable to calculations 
based on any symmetry restoration or generator coordinate but is limited to EDFs depending only 
on integer powers of the normal and anomalous density matrices. Eventually, the method is formally 
illustrated for particle number restoration and is specified to configuration mixing calculations based 
on Slater determinants. 

PACS numbers: 21.10.Rc, 21.60.Ev, 71.15.Mb 

Keywords: Energy density functional, configuration mixing, symmetry restoration, self-interaction, self- 
pairing 



I. INTRODUCTION 

The nuclear Energy Density Functional (EDF) ap- 
proach is the microscopic tool of choice to study medium- 
mass and heavy nuclei in a systematic manner The 
EDF approach used in nuclear physics has formal simi- 
larities with Density Functional Theory (DFT) [3,[l,[l,i, 
H, S 111 which provides a formal framework to obtain the 
exact ground-state energy and one-body density of elec- 
tronic many-body systems in condensed-matter physics 
and quantum chemistry . However, and even if it is of- 
ten referred to as nuclear DFT 0, US, [H, [H, Q [H , the 
nuclear EDF method has deeply rooted conceptual dif- 
ferences with standard Density Functional Theory; e.g. 
Ref. [3. 

One of the major challenges of the nuclear many-body 
problem is that atomic nuclei are self-bound strongly cor- 
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related composite systems that tend to display strong 
collective modes and individual excitations on the same 
energy scale. This makes the modeling of nuclei rather 
different from that of electronic systems in external po- 
tentials to which DFT is heavily applied. The underlying 
assumption to nuclear EDF approaches is that correla- 
tions can be divided into different classes that can be in- 
corporated through successive refinements of the method. 
While short range in-medium correlations are subsumed 
into a suitable energy density functional, long-range cor- 
relations associated with collective modes must be in- 
corporated more explicitly. On this basis, two different 
levels of EDF calculations coexist in nuclear physics. 

On the first level, inappropriately called self-consistent 
mean-field theory, the EDF is constructed from a one- 
body density matrix that corresponds to a single prod- 
uct state. This reference state is either taken as a Slater 
determinant, abusively referred to as Hartree-Fock (HE) 
state, or as a quasi-particle product state, abusively re- 
ferred to as Hartree-Fock-Bogoliubov (HFB) state. In the 
present work, this level of description will be referred to 
as the Single Reference Energy Density Functional (SR- 
EDF) method. Within the SR-EDF approach, static cor- 
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relations associated with collective modes are incorpo- 
rated through the use of a symmetry-breaking product 
state. 

On the second level, several of those symmetry- 
breaking product states are mixed in the framework of 
the Generator Coordinate Method (GCM) with the aim 
of restoring broken symmetries, allowing for fluctuations 
of the order parameters of broken symmetries, or both. 
Such a Multi-Reference (MR) EDF approach allows one 
to incorporate correlations associated with large am- 
plitude collective motion beyond the static correlations 
which are easily treated within the SR-EDF method. 
The energy functional at play in MR-EDF calculations 
depends on the transition density matrices constructed 
from all possible pairs of product states that enter the 
MR set. 

Until now, there exists neither a rigorous formal frame- 
work underlying MR-EDF calculations nor a guidance 
from DFT. Indeed, the nuclear MR-EDF method is 
outside the scope of existing implementations of the 
Hohcnberg-Kolm theorem, even of those dealin g sp ecif- 
ically with symmetry breaking issues 
As far as motivating MR-EDF calculations based on the 
GCM, hybrid approaches used in electronic systems that 
combine a density functional for diagonal energy ker- 
nels and (sometimes scaled) matrix elements of the bare 
Coulomb force for off-diagonal ones [2l| have to be ruled 
out because of the different nature of the nuclear force 
and the more complex correlations it induces in the nu- 
clear medium. There is a consensus among practition- 
ers that the SR-EDF should be recovered for any diag- 
onal energy kernel appearing in the MR-EDF and that 
the non-diagonal kernels should be extrapolated from the 
SR-EDF on the basis of the Generalized Wick Theorem 
(GWT) [2^ which provides an efficient tool to evaluate 
the matrix elements of any operator between two differ- 
ent product states. For diverse reasons, however, SR and 
MR-EDFs do not reduce to matrix elements of a Hamil- 
tonian operator between many-body states. As a con- 
sequence, the procedure outlined above to generate the 
MR-EDF is not free from problems and inconsistencies. 

The present article is the first paper of a trilogy were 
we address one of the problems arising from the current 
lack of a stringent and consistent constructive frame- 
work for the MR-EDF method, namely the appearance 
of divergences and finite steps that have been recently 
identified in the context of Particle Number Restoration 
(PNR) [H, 113, m, m, [13. We propose here a general 
cure for any type of configuration mixing. In Ref. psj 
the method is applied to PNR while in Ref. [2^ the spe- 
cific case of an EDF depending of non-integer powers of 
the density matrix is discussed. The spurious steps and 
divergences arise typically as one scans the symmetry- 
restored energy as a function of a certain degree of free- 
dom; e.g. looking at the PNR potential energy surface of 
a nucleus as a function of axial quadrupolc deformation. 
For reasons that will become clear below, pure PNR is 
the only case where the divergent contributions to the cn- 
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ergy functional can be traced back rather easily although 
no practical solution to this problem has been proposed 
so far. It is important to realize, however, that similar 
problems may appear for any MR calculation performed 
within the EDF framework, although they will be hidden 
in the formal expressions and may be masked by the nu- 
merical representation. A pathology that resembles the 
one identified for PNR has recently been seen in EDF cal- 
culations, without explicit treatment of pairing, aiming 
at restoring angular momentum (30j . whereas a similar 
problem was identified much earlier in the GCM-type 
mixing of zero- and two-quasi-particle states (3lj | . 

It is essential to state that these difficulties are inherent 
to the EDF formulation of MR calculations and do not 
exist when the energy relates to the strict average value 
of a Hamiltonian in a wave-function [l^, US]- As 
the major part of the paper consists of lengthy formal 
manipulations that will be easier to follow if anticipating 
the results, we now summarize the three goals of the 
present work: 

(i) We want to introduce a formal framework that al- 
lows the unambiguous identification of terms in the 
EDF that are responsible for divergences and finite 
steps in MR-EDF calculations, irrespective of the 
nature of the configuration mixing involved. As 
will be shown below, finding the canonical basis 
of the Bogoliubov transformation which links two 
HFB states allows us to do so. In that basis, the 
non-diagonal energy kernel associated with those 
two vacua and which enters the MR energy can 
be constructed on the basis of the Standard Wick 
theorem instead of the Generalized Wick Theorem. 
This will be the key to our analysis. 

(ii) Guided by the strict Hamiltonian case, several 
sources of problems including those responsible for 
divergences and steps are unambiguously identified 
within the EDF context. 

(iii) We provide a minimum but general solution to the 
problem of divergences and steps in MR-EDF cal- 
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culations for a particular class of energy function- 
al. To that end, we propose a correction that 
removes exactly the difference between the GWT- 
based and the SWT-based non-diagonal energy ker- 
nels. The latter correction preserves the continu- 
ity condition between SR and MR-EDFs. Finally, 
we provide explicit expressions for particle-number 
restoration and mixing of Slater determinants. 

The paper is organized as follows: Section [ill briefly 
reviews general aspects of nuclear EDF approaches. Sec- 
tion |llT] is devoted to the actual introduction to SR and 
MR-EDF methods. First, the relevant similarities and 
key differences between Hamiltonian-based calculations 
and EDF ones are worked out. Second, the connection 
between the SR functional and the non-diagonal energy 
kernels at play in the MR formalism is discussed. In 
Section IIVI we demonstrate that the usual definition of 
non-diagonal energy kernels based on the GWT leads to 
ill-defined calculations when performed within the EDF 
approach. In Section|V]we propose a correction removing 
the spurious contributions to the energy functional that 
is applicable to any MR-EDF calculation. The correction 
is discussed in detail in Sections lVII and IVlIl for two appli- 
cations of current interest. Conclusions and perspectives 
arc given in Section IVlIII 



II. GENERAL ASPECTS OF THE NUCLEAR 
EDF APPROACH 

The nuclear energy density functional integrates out 
short-range correlations associated with non-collective 
excitations on a given vacuum, in particular the strong 
tensor correlations from the vacuum nucleon-nucleon 
force, thereby introducing higher-order density depen- 
dencies into the energy functional. The diagrams that 
are summed up in the particle-hole and particle-particle 
channels arc different [s^, which strongly suggests 
the use of different effective vertices in those two chan- 
nels of the EDF. Lacking an explicit connection to first 
principles, existing functionals are constructed using the 
symmetries of the nuclear Hamiltonian as a guiding prin- 
ciple whereas their parametcrizations are adjusted phe- 
nomenologically [ij. 

Allowing the reference state to break the symmetries of 
the eigenstates of the underlying Hamiltonian is a way to 
incorporate static long-range correlations associated with 
collective modes, as for example deformation and pair- 
ing, with very moderate effort. However, the breaking 
of symmetries (translational, rotational, parity, particle 
number, to name the most current ones) forbids a triv- 
ial connection of the nuclear SR-EDF formalism to the 
Hohenberg-Kohn (HK) theorem which is the foundation 
of DFT [l^ . Indeed, the density that minimizes the exact 
HK energy functional must reflect the symmetries of the 
exact ground state of the system. In fact, the appearance 
of symmetry-breaking solutions in nuclear EDF calcula- 
tions underlines two important elements (i) it is crucial 



and numerically not too difficult to grasp the most impor- 
tant static correlations using rather simple approximate 
functionals and a single-determinantal reference state, at 
the price of violating the HK theorem (ii) kinematical 
correlations associated with the corresponding symmetry 
modes (Goldstonc modes) as well as the correlations due 
to the ffuctuation of their order parameters are extremely 
difficult to incorporate into a single-determinantal ap- 
proach. In other words, correlations associated with 
highly non-local processes such as large-amplitude col- 
lective motions can hardly be described within a SR ap- 
proach based on a standard, nearly-local EDF. 

It is possible to improve on the SR level by extend- 
ing the EDF formalism to a Multi-Reference frame- 
work, inappropriately referred to as beyond-mean-field 
calculations in the literature. In a MR calculation, a 
set of (usually) non-orthogonal product states of ref- 
erence are combined to construct a more general en- 
ergy functional which involves the computation of non- 
diagonal norm overlaps and energy kernels between all 
pairs of vacua. Such MR-EDF calculations are inspired 
by projection techniques and the Generator Coordinate 
Method (GCM) whose formalisms are well documented 
for Hamiltonian-based calculations [H, [H, [s^, [13 . Sym- 
metry restoration and GCM provide an efficient way to 
incorporate both collective and single-particle dynam- 
ics into a coherent quantum-mechanical formulation. In 
particular, the GCM allows the inclusion of correlations 
associated with the fluctuation of the order parameters 
of the broken symmetries. Group theoretical consider- 
ations [1^, the Random-Phase- Approximation limit of 
GCM-type EDF calculations [s^, as well as the require- 
ment of continuity between SR and MR functionals dis- 
cussed in the introduction can be used to constrain the 
form of the energy kernels to be used in a MR formal- 
ism. In particular, it can be shown that the latter energy 
kernels have to be constructed exclusively from transi- 
tion normal and anomalous one-body density matrices, 
cf. Eqns. (|14m6p below, between the two reference states 
involved. These transition density matrices, however, are 
at the origin of a major difficulty that arises in existing 
implementations of MR EDF calculations. 

It is noteworthy that the description of particular sys- 
tems and phenomena in electronic systems also calls for 
extensions of the Kohn-Sham implementation of HK- 
DFT to deal with symmetry-breaking solutions [13, [H, 
[igl . [20j . The latter extensions are still connected to the 
HK theorem and differ from MR-EDF calculations per- 
formed in nuclear physics to restore broken symmetries. 



HI. FROM SINGLE-REFERENCE TO 
MULTI-REFERENCE FUNCTIONAL 
FRAMEWORKS 

We consider the nuclear Hamiltonian under the form 
^ = Y1 + ^ XI ^ii^i ^-t o-^ «i , (1) 

ij ijkl 
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which, for the sake of transparency of the chain of ar- 
guments below, has been limited to the sum of the ki- 
netic energy and a two-body interaction ui2, represented 
through its antisymmetrized matrix elements Vijki- For 
the moment we also do not distinguish between protons 
and neutrons. The extension to the case of two nu- 
clcon species and three-body or even multi-body forces 
of higher order is straightforward (although cumbersome) 
and will not influence our conclusions. 

It is important to stress at this point that, in the 
present paper, H defined above strictly refers to a "vac- 
uum" nuclear Hamiltonian operator; e.g. its two-body 
part reproduces nucleon-nucleon scattering data up to 
relevant laboratory energies. Later in this section we re- 
fer to effective Hamiltonians, either in the context of a 
truncated single-particle space or because they re-sum 
correlations through an explicit dependence on the den- 
sity of the system. It happens that replacing the vacuum 
Hamiltonian H hy a, genuine operator defined within a 
truncated single-particle space would not modify the fol- 
lowing discussion regarding the technical issue addressed 
in the present paper, whereas using a density-dependent 
"Hamiltonian" is at the origin of some of the problems 
we are concerned with. 



A. Single-reference framework 

In the SR approach, the reference state is taken to 
be an independent-quasi-particle state 1^*0)1 most conve- 
niently written as a quasi-particle vacuum 



through 



(2) 



where Co is a complex normalization coefficient. The 
set of quasi-particle operators {a^, a+} is obtained from 
the complete basis of single-particle operators {ai,af} 
through a Bogoliubov transformation characterized by 
the matrices {U^,V°) 



(3) 



The index "0" appearing in |$o) as well as in J/" and 
V° denotes an ensemble of collective coordinates that 
fully characterizes the product state |<&o), either regard- 
ing broken symmetries or the average value of a con- 
straining operator. This label is of no explicit impor- 
tance at the SR level but will become mandatory for MR 
calculations. 



1. Strict HFB approach 

In the strict HFB approach [40,], the reference state 
|$o) approximates the actual ground-state of the system 
in such a way that the approximate energy is obtained 



Eo 



(^o|g|^o} 
(*o|*o) 



(4) 



Using the Standard Wick theorem (SWT) [4l|, Eq. (U) 
can be expressed as 



^00 „oo 
ijkl Pki Pi J 



ijkl 



+iE%'^'«°r'^°?, (5) 



ijkl 

where p*^" is the intrinsic (SR) normal one-body density 
matrix and k*^" is the intrinsic anomalous density matrix 
(also called pairing tensor) 



00 ^ (^0 1 1^0) 

~ ($o|$o) 

no ^ (^o|aja»|^o) 
«y - 



(6) 
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OttO+\ 



(8) 



associated with the reference state |<&o)- The normal den- 
sity matrix is hermitian p^j* = /5^° * , whereas the anoma- 



^00 



lous density matrix is skew symmetric K^j = 

The energy Eq given by Eq. ([5]) can be seen as the 
particular energy functional £^ [p™, k*^" *] that re- 
sults from taking the expectation value of the Hamil- 
tonian in the product state. This is what is referred to 
as the strict HFB approximation. The minimization of 
E [/9°°, with respect to all independent degrees 

of freedom (p^'° as weh as pff*, K^f, and K^f* for 
j < i), under the constraints that the reference state 
|$o) remains a quasi-particle vacuum and that the av- 
erage particle number ($o|-/V|<I>o)/($o|^o) has the fixed 
value TVo, leads to the HFB equation that determines 
{U°,V°) and the densities. 

So far, we have not discussed the characteristics of 
the nuclear Hamiltonian H . Traditional nucleon-nucleon 
hard-core potentials are not perturbative [i^l, which 
makes any attempt to use the strict HFB approximation 
to the exact nuclear many-body problem useless. The re- 
cently proposed soft-core interactions on the other hand 
seems to make the nuclear many-body problem pertur- 
bative [i^. Still, one must go beyond lowest order to 
obtain close to converged results and the strict HFB ap- 
proximation is not quantitatively viable. The situation is 
further complicated by the necessity to treat three-body 
and maybe higher-body forces in the nuclear many-body 
problem [is, 44 , 45), Ho, [13 ■ As a consequence, one often 
resorts to an effective Hamiltonian that implicitly incor- 
porates correlations brought by the physics outside the 
model space used in a given calculation. Traditionally, 
there are two major classes of effective interactions for 
self-consistent mean-field calculations that represent two 
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different philosophies and strategies to reduce the nu- 
clear many-body problem to a tractable number of rele- 
vant degrees of freedom. There are (i) approaches that 
can be seen as lowest-order approximations to the inter- 
acting shell model [i^, formulated in terms of effective 
(often schematic) Hamiltonians in a strict HFB approach 
using a schematic shell-model space of a very few spher- 
ical harmonic-oscillator j shells [s^ (ii) approaches that 
use the full model space of occupied particles together 
with a density-dependent effective interaction, which are 
what is nowadays recognized as an approximation to a 
more general SR-EDF formalism. 

2. EDF approach 

The EDF approach builds upon the last comment. In 
the context of an EDF formalism, the reference prod- 
uct state |$o} used to construct the one-body normal 
and anomalous density matrices (OS]) is to be seen as 
an auxiliary state which is not meant to be a good ap- 
proximation of the actual ground-state wave function, 
in a similar manner as in DFT for electronic systems 
i, H, i, 0, i, B m • The EDF and the HFB-like equation 
which results from its minimization contain higher-order 
correlations than those provided by the strict Hartree, 
Fock and Bogoliubov diagrams written in terms of the 
vacuum interaction. It is important to note that, as op- 
posed to the strict HFB method, the SR-EDF approach 
does not rely on the Ritz variational principle in the sense 
that the minimal energy obtained from an approximate 
energy functional may be below the actual ground-state 
energy of the system. 

Thus, instead of the expectation value of H, the 
starting point of the method is an energy functional 
£ of the normal and anomalous density 

matrices. Considering the simple case of a bilinear func- 
tional in order to stay formally close to the strict HFB 
approach with a two-body interaction, a typical EDF can 
be written as 

£:[p™,K°°,K™*] = SP + £PP + £'"' 

Pji + 2 ^'ijkl Pki Pi] 
ij ijkl 

+ iE^SX.'^??*-2?(9) 

ijkl 

In Eq. the matrix elements of the effective vertex v^p 
might or might not be antisymmetric; e.g. due to the 
Slater approximation used to treat the exchange contri- 
bution from Coulomb or by dropping specific terms of 
the Skyrmc functional Due to the antisymmetry of 
K, however, only the antisymmetrized part of the vertex 
is probed in the last term of Eq. ^ and one can always 
take v'^'^ to be antisymmetric, which we will do here. In 
any case, and even though Eqs. ^ and ^ look very sim- 
ilar, yPP and v'^'^ should not be seen as interactions in the 
real sense and arc likely to be different. Of course, they 



are in principle related to the original vacuum interaction 
but, in the absence of a constructive framework, the link 
remains implicit. Popular energy density functionals for 
calculations along these lines [l|] are the non-relativistic 
Fayans [l2|jGogny [sO] and Skyrme [5l[ as well as the 
relativistic [53| ones. 



B. Multi-reference extension 

Because the nuclear Hamiltonian is invariant un- 
der certain symmetry transformations, the independent 
quasi-particle state |$o) should be an eigenstate of ap- 
propriate linear combinations of the infinitesimal gener- 
ators. However, the simultaneous preservation of sym- 
metries and the inclusion of long-range correlations asso- 
ciated with large-amplitude collective motions is impos- 
sible within the strict HFB formalism and remains un- 
tractable within a SR-EDF formalism. As a result, the 
product state |$o) is permitted to spontaneously break 
the symmetries of the true eigenstates to lead to energet- 
ically favored solutions; an ambiguity of mean-field and 
EDF approaches for which Lowdin coined the notion of 
" symmetry dilemma" in Ref . (53| . 

To solve that ambiguity, one resorts to a multi- 
reference extension of the method [ij . The starting point 
consists of a finite set of product states {I'I'o); G MR}, 
where {0 S MR} denotes different realizations of the 
set of collective coordinates that characterize a reference 
state. If focusing on symmetry restoration, the states be- 
longing to the MR set correspond to one another by the 
application of a transformation of the symmetry group 
under consideration. However, the label of the product 
states in the MR set can also refer to different values of 
the order parameter of the broken symmetry; the goal 
of mixing such states being to incorporate correlations 
associated with fluctuations of that order parameter. In 
such a case, the operator that links two vacua in the set 
is not known analytically. 



1. Strict projected-GCM approach 

In the strict projection-GCM approach, one constructs 
the projected-GCM state from the product states in the 
MR set through 

I*')- E /ol*o), (10) 

{0}eMR 

where the determination of the weight functions /q is 
discussed below. The index k denotes that not only the 
ground state, but also excited states can be extracted, 
either by projection, by diagonalization, or both. In the 
strict projection-GCM approximation, the many-body 
energy is approximated by the average value of the nu- 
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clear Hamiltonian in the projcctcd-GCM wave-function 
, _ (vl/^|i/|vl/fc) 



(11) 



(12) 



where the set of energy i?[0, 1] = (^o\H\^i) / {^o\^i) and 
norm ($o|^i) kernels constitute the basic ingredients of 
the method. Note that Eq. ([T^ is nothing but the en- 
ergy in the (Hamiltonian-based) Generator Coordinate 
Method (GCM) [H, [H. From a formal point of view, 
symmetry restoration is a special case of the GCM. In- 
deed, in that case j^'*^) can be expressed as a projector 
acting on one of the symmetry-breaking product state in 
the MR set. Whereas in the GCM the weight functions 
f'' arc determined variationally from the Hill-Wheclcr- 
Griflin equation [H^, \EM j they are dictated by properties 
of the relevant symmetry group when performing sym- 
metry restoration^. 

The energy kernel E[Q, 1] can be evaluated through the 
use of the Generahzed Wick Theorem (GWT) ^ 



^01 ^01 
"ijkl Pki Plj 



ij ijkl 



ijkl 



It is essential for the present work to note that Eq. ([Te 
is exactly of the same functional form as Eq. ([5]), except 
that intrinsic density matrices {p'^*', k*^", k'^^*} have been 
replaced by transition density matrices defined as 



($o|a+a,|$i) 

($o|*i) 
($o|aja,|$i) 

($o|*i> 
($o|a+a+|$i 
($o|*i) 



(14) 
(15) 
(16) 



where k^- * is in general not the complex conjugate of k] 



anymore. Given Eq. p^ . it is clear that any "diagonal" 
energy kernel, i.e. obtained for |$i) = |<i>o)j is equal to 
the strict HFB energy, Eq. ([5]), so that the continuity re- 
quirement between the SR energy and MR energy kernels 
stated in the introduction is trivially fulfilled. 

As in the strict HFB approximation, a strict projected- 
GCM calculation performed in terms of the vacuum nu- 
clear Hamiltonian does not lead to quantitatively satis- 
factory results for all observables of interest. Still, the 
restoration of the Galilean invariance of spherical HF 



states |56l.l57j] was studied within such a scheme; i.e. omit- 
ting the density-dependent part of the Gogny interaction 
and keeping exactly all exchange terms in the computa- 
tion of MR kernels. In fact, MR calculations employ- 
ing Hamiltonian operators are rather performed within a 
limited valence space, in the spirit of the interacting shell 



6l| or 



model, cither of projected-GCM type [58|, l59|, l6C 
using more elaborate schemes to select the set of reference 
states like in the MONSTER and VAMPIR approaches 

MMM- 



2. EDF approach 

To perform quantitatively relevant calculations using 
the full model space of single-particle states, one has to 
turn to an energy density functional variant of MR cal- 
culations. Guided by Eq. the many-body energy is 
defined in the MR-EDF approach by 



^,^E{o,i}emr/o*/i'^[0,1] (*o|<I>i) 



X]{0,l}eMR, /o 



./iMl'ol^i 



(17) 



which now depends only implicitly on the projected- 
GCM state of Eq. dTO]). Indeed, the energy £^ is not the 
average value of or any genuine operator^, in \^^) 
and it cannot be re-expressed into Eq. (fTT|) . Rather, £^ 
is to be seen as a more general functional of all product 
states belonging to the MR set {|<I>o);0 G MR] as each 
energy kernel 5[0, 1] is a functional of a pair of states 
{|$o);|*i)} in the set. 

Although Eq. p?)) provides the general ground to cal- 
culate the energy of the system within the MR-EDF for- 
malism, there remains the question of how the energy 
kernel £[0, 1] should be constructed and what its connec- 
tion to the SR energy functional £[p'^^ , k^^ , *] intro- 
duced in the previous section is. Guided by the structure 
of the expressions in the strict HFB and GCM theories, 
Eqs. (|5l) an d (flSt, p ractitioners of nuclear EDF methods 
[25l . |39[ m, l66l. l67l. [68j have used the natural extension 



^[0, 1] = 



01 ^01 ^10*1 



to define the MR energy kernel 



entering Eq. (|17p from the SR energy functional. Using 
the bilinear SR functional of Eq. ([9]) for illustration, this 
leads to a kernel of the form 



GWT 



[0,1] 



ijkl 



IV^-KK ^.01*^.10 /ION 

4 2-^ "ijkl '^ij '^kl , \^°) 



ijkl 



which is labeled as "GWT" since the definition £\0, 1] 



01 ^01 ^10* 



] amounts to using the generalized Wick 



^ A density-dependent operator, that is an operator that depends 
^ The weight functions are fully determined only if the symmetry on the solution, is not considered in the present work as a genuine 

group is abelian. operator. 
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theorem as a motivation to define the MR energy kernel 
from the SR-EDF. 

Calculations along these lines using realistic EDFs 
have been performed for the restoration of particle num- 
bers (25l . [68j and ang ular momentum of axially-deformed 

of angular momentum o ' 



HFB states 1691, |7i 
cranked HF states 3^^ 



■ 73ll. and of parity of octupole de- 
formed HFB states Some of those calculations 
also included a GCM-type configuration mixing along a 
collective degree of freedom related to spatial deforma- 
tion. 

Other prescriptions than the simple replacement of 



pTi with 



Pit 



to unreaUstic results [39[. The prescription £[0,1] = 
£[p^^ , , k}^ *] ensures that the continuity requirement 
between the SR and MR levels of description is fulfilled. 
For symmetry restorations, it can also be shown that 
the functional should only depend on transition densi- 
ties [H. 

However, the GWT-based prescription to define f [0, 1] 
causes seriousproblems as was recently illustrated for 
PNR in Refs. l25l.[26l | and confirmed from a different point 
of view in Ref. |28l|Tafter earlier warnings [H, [13, [3l| were 
not recognized. The most spectacular manifestation re- 
lates to the appearance of divergences in the PNR energy 
when two vacua in the MR set are orthogonal, which in 
pure PNR applications happens for a relative gauge angle 
of 7r/2 when a pair of sing le-p article states has the occu- 
pation probability of 1/2 [2^. While this does not 
any problem in the strict PNP-HFB method 
this makes PNR to be ill-defined within the MR-EDF 
approach. In fact, the energy functional does not only 
contain divergences when encountering orthogonal vacua, 
but also finite spurious contributions even when the two 
vacua are not orthogonal [2^ as will be illustrated later 
in the present work. It is the goal of the following sections 
to introduce a method that allows the identification and 
removal of these spurious terms from any type of MR- 
EDF calculations. 

It is noteworthy that, to the best of our knowledge, 
the continuity requirement that is always used as a guid- 
ing principle to construct nuclear MR energy kernels 
is never enforced in MR calculations for electronic sys- 
tems. There, hybrid approaches that use a density func- 
tional for diagonal kernels and the (sometimes scaled) 
bare Coulomb force for non-diagonal kernels are used in- 
stead [l^ [10, [t^, [tJ , an approach that is facilitated by 
the fact that the correlation and exchange parts of the 
energy are a mere correction to the direct Coulomb in- 
teraction in most Coulomb systems. 



countered in the latter. As a result, we first concentrate 
on the strict projected-GCM method and compare the 
computation of the energy kernel {'^o\H\(^i) / {'^o\<^i) as 
obtained from both the GWT and the SWT in a suitable 
basis. 



A. Notation and preliminary discussion 

In this section we summarize the elements which will 
be needed below for the analysis of MR-EDF calcula- 
tions. The key ingredient will be to find a (quasi-particle) 
basis that allows the computation of the energy kernel 
($o|i^|$i>/('J'o|$i> in terms of the SWT instead of the 
commonly used GWT. The question if and how this is 
possible has been addressed already from various per- 
spectives in the past [3(3, 58, ,78, . 

1. We write the two vacua involved as products of the 
two corresponding (different) sets of quasi-particle 
operators, denoted by and (3^^ respectively. 



i*o) = Co n^-io)' 



(19) 



The moduli of the two complex constants Co and Ci 
are fixed by the normalization of the states, while 
their phases can be set to any value. 

2. We introduce two sets of matrices (C/",!^") and 
{U^ ,V^) that represent the Bogoliubov transfor- 
mations between an arbitrary, but common, single- 
particle basis {02,0^} and the two sets of quasi- 
particle states 



a 



i 



(20) 
(21) 



3. One realizes that the transformation that expresses 
the quasi-particle operators associated 
with in terms of the quasi-particle operators 
{a^,Q;+} associated with [$o) is canonical, which 
means it has the form of a general Bogoliubov 
transformation 



(22) 



IV. ENERGY KERNELS IN MR APPROACHES 

Although the main point of the present work is to stress 
the differences between strict HFB/projccted-GCM ap- 
proaches and EDF ones, the analysis of the former does 
offer the key to the understanding of the problems en- 



with the matrices A and B satisfying (see Ap- 
pendix E of Ref. [3^) 



A EE [/0"f/i + \/0"F\ 
B = V°^U^ + U°^V^ . 



(23) 
(24) 
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To avoid any misunderstanding in what follows, we 
stress that the Bogoliubov transformation rep- 
resented by the A and B matrices has nothing to 
do with pairing correlations. As a consequence, its 
ingredients have a particular meaning and inter- 
pretation that might seem unusual and unexpected 
(e.g. see Sec. lVII|) . 

We also stress that no restriction or symmetry is 
presently imposed on the reference states {|<i>o); € 
MR}. They might break time-reversal invariance 
and have an odd number-parity (sol . [Slj (describing 
odd nuclei), be two-quasi-particle states (describ- 
ing diabatic excitations), or higher quasi-particle 
states. The only restriction is that all reference 
states belonging to the MR set have the same def- 
inite number parity. 

4. The matrices A and B defined through Eqs. (|22j|24p 
can be decomposed into a sequence of three sim- 
pler transformations thanks to the Bloch-Mcssiah- 
Zumino (BMZ) theorem [sS, [H] 

A = DAC, B = D*BC. (25) 

where D (C) is a unitary transformation which 
only mixes quasi-particle creation/annihilation op- 
erators {a+}/{Q:,y} ({/3^}/{/3p}) among each other, 
whereas A and B represent a special Bogoliubov 
transformation. The transformations D and C 
of the decomposition (|25p introduces two interme- 
diate quasi-particle bases {a^,a'^} and {f3p,,fi'^} 
with different properties from the two original ones 
{a,y, a'^} and {(3^,, As a quasi-particle vacuum 
is invariant under a unitary transformation among 
quasi-particle creation/annihilation operators that 
define it, one could initially choose |$o) and |$i) 
such that C = D. In practice, however, this is usu- 
ally not the case as |$o} and |$i) are constructed 
independently from each other and have to be seen 
as given inputs to a MR calculation. Finding the 
transformation that allowed to set C = D would 
constitute an effort similar to performing the de- 
composition (ps)) . 

The advantage of the Bloch-Mcssiah-Zumino de- 
composition of the transformation (A, B) becomes 
clear when looking at the two intermediate sets of 
quasi-particle operators a+ and /3+ 

a+^Y. ' /5+ ^ E ■ (26) 

The two product states |$o) and |$i) arc still vacua 
for the two sets of quasi-particle operators {olh, aj^} 
and {f3y,(3^}, respectively. In addition, the trans- 
formation expressing the latter in terms of the for- 
mer takes the simple BCS-like form of a special 
Bogoliubov transformation 

/3+ = A^^ 5+ + 5p . (27) 



The matrices A and B are 2x2 block-diagonal 
(after a suitable rearrangement of the sequence of 
indices), with the non- vanishing blocks A(p) and 
B{p) being defined as 

(28) 

where App = App and Bpp ~ —Bpp. The block- 
diagonal structure of A and B shows that the quasi- 
particles {dt^^a'^} come in conjugated pairs {v.,v), 
whose nature will be discussed below. As usual, 
the block structure can be used to formally sepa- 
rate the corresponding quasi-particle basis into two 
halves, even though it is not necessary at this point 
to specify which quantum number is used to distin- 
guish them. For the rest of the paper, we introduce 
the convention that the label p refers to > while 
p refers to < 0. Thus, when v ~ p{p), then 
D=p{p). 

5. Standard BCS techniques allow us to express |$i) 
in terms of |3>o) through 

= Coi n i^*PP + Kp^t ^t) l*o) , (29) 

where Cqi is a complex normalization factor whose 
properties will be discussed in Sec. IIVBI 

There are a few further remarks to be made before we 
continue. When A* ^ for allp, Eq. (^5)) is nothing but 
Eq. (E.39) of Ref. ^al- However, the form above is also 
valid when |$o) and j^i) are orthogonal, i.e. when at 
least one pair (p,p) is such that A*p = (i.e. B*^ = 1). 
This situation corresponds to Eq. (E.38) of Ref. with 
n even and "blocked" quasi-particles coming in pairs'^ 
{p,p). One such case occurs when |$i} is obtained as a 
two-quasiparticlc excitation {po,Po) on top of |$o)- This 
amounts to having {App ~ l,Bpp = 0) for all pairs but 
the pair {pq,Po) for which {Ap„pg ~ 0, Bp^p^ = 1). In fact, 
there are two different situations to consider depending 
on the properties of the matrix D. If the latter is the unit 
matrix, the two-quasiparticle states created are conju- 
gated in the original basis {a^,a~^} used to characterize 

$o)- However, other two-quasi-particle excitations on 
top of |$o) can be generated by considering a non-trivial 
D transformation. 

The above procedure to express |$i} with respect to 

$o} is close to the method discussed in Ref. [7^, al- 
though the latter relies on the Thouless representation of 
quasi-particle vacua [s^ l which requires to take a more 
explicit care of orthogonal states. Refs. [H, [tI, [7§] are 



^ The conjugated pairs arc not necessarily related 

through time-reversal transformation. 
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also dedicated to finding a canonical representation ap- 
propriate to non-diagonal overlaps. However, they focus 
on the more difficult task of finding a basis where both 
vacua are simultaneously put into a canonical form. In 
the present work, we are only interested in putting the 
Bogoliubov transformation {A, B) describing the transi- 
tion between the two vacua into a canonical form but not 
the vacua themselves. 



B. Overlaps and one/two-body contractions 

The form given by Eq. (|29[) is very convenient to ex- 
press the overlap and the matrix element of operators 
between the states |$o) and |$i}. First, let us precise 
the nature of the normalization constant Cqi. As already 
said, the two vacua are usually taken to be normalized 
through an appropriate choice of |Co| and |Ci| in Eq. 
Calculating 

($i|$i) = |Coi|' ($o|$o) n (l^Pfl' + l^wl') (30) 

p>0 

from Eq. (HH) and using that |^ppp -I- |-BppP = 1 for aU 
p, the normalization of the two vacua leads to |Coip = 
1. However, the latter result does not fix the phase of 
Coi which remains when calculating the norm overlap 
between the two vacua 

($o|*i) = Coi n Kp = Coiy^dct(A*). (31) 

Eq. ((3T|) differs from the well-known Onishi-Yoshida for- 
mula [sHi by that phase factor Cqi which depends on the 
chosen convention, that is, on the explicit form of the 
states |$o) and |$i) taken as an input of the MR cal- 
culation. We will show exphcitly in Sec. |Vl] how this 
works for PNR where the phase Cqi can be obtained 
analytically for all pairs of vacua involved in the cal- 
culation. Unfortunately, this is not true for more in- 
volved MR calculations, as for example when superpos- 
ing particle-number and angular-momentum restorations 
together with GCM-type configuration mixings. In such 
a case, one must desi gn n umerical methods to follow the 
phase Coi [H, [zl, HlsJ . 

For reasons that will become clear in a moment, we 
define for p ^ q the quantities 

($o|$i,p) EEE Coi n ^pV' (32) 

p'>0 

($0 1*1,^,9) = Coi n ^^'P" (33) 

p'>0 

using the notations of Ref. [H. In Eqs. p2ll33| . 
(<i>o|<i>i,p) and ($o|^1jP, ?) remain unchanged when sub- 
stituting p and/or qhy p and/or q in such a way that one 
can refer to ($o|$i,i^) and ($o|^i, m)- Note that the 



latter overlap does not need to be defined iov v = ji or 
V — as. it cannot appear in such cases. Indeed, this 
would correspond to removing a conjugated pair twice 
from |$i), which would lead to 

\<^i,v,v) oc Q;,.ap|$i,j/) . (34) 

If needed, such particular cases can be trivially taken into 
account by extending the notation introduced in Eq. (|33|) 
through 

($ol*i,i^,J^) = (*ol*i,J^,'>) =0 . (35) 

Starting from Eq. it is now straightforward to cal- 

culate one-body contractions in the quasi-particle basis 
{a^, a+} using the SWT. Using the notations introduced 
in Eqs. [3T]and[32l one finds that 

($o|fij5p|*i) = ($o|5+ci+|$i) =0, 

($o|5.5p|$i) = <5p^b;,($o|$i,i^)- (36) 

All two-body contractions are obtained in the same man- 
ner and can be written in a compact way using the no- 
tation introduced in Eqs. [551 and [55l e.g. 

= (5p^(5^a-B;^B;^($o|*i,i^,7) 
5p.s B*,^B*,^ ($o|$i,i^, m) 
+<5,5<5p^5*,5*^($o|$i,t^,M)- (37) 

C. More convenient bases 

In order to take advantage of the results obtained in 
the previous section, one has to express the single-particle 
creation and annihilation operators {ai,af} in terms of 
the quasi-particle operators {Q:i/,a+}. Introducing the 
matrices = U'^D and V'^ = V^D, we have 

Many formulas derived below will simplify by introduc- 
ing the quasi-particle wave-function associated with the 
quasi-particle operators {5,^,5+} 

(!::;)■ 

where the choice of labeling the lower component \ipp) 
with the quantum number of the conjugated state in the 
other half of the basis is dictated by convenience, as will 
become clearer in the PNR case discussed in Sec. IVII 

The upper and lower components can be expressed 
in terms of the arbitrary single-particle basis {ai,af} 
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through 



i 



(40) 
(41) 



which are at variance with those expressed in terms of ro- 
man indices [i, j,k,l) which relate to the initial arbitrary 
basis {ai,af}. 



We recall that C/° and V° introduced above involve a uni- 
tary transformation D determined by the Bloch-Messiah- 
Zumino decomposition of the Bogoliubov transformation 
{A,B). These definitions will allow us to use matrix el- 
ements of the effective vertices v'''' and v'^'^ written in 
the mixed basis of upper and lower components of the 
quasi-particlc states {|(?I>i/), Iv'p)}, e.g. 



pp 



EyO T>0 fj-a frO -PP 

ijkl 



(42) 



D. Energy kernel from the SWT 



Thanks to Eqs. ([36]) and ((37|) . the contribution 
{'^o\vi2\'^i) / to the energy kernel £;svft[0, 1] cal- 
culated through the application of the SWT reads as 



i?g;,^[o, 1] + e^wtIq, 1] ^\ E + J E 



Ufl \ U I / \ I I / 

lyfi \ I / \ I / 

1 V 71 B* B* ($o|$i,i/,m) 1 V- f5* (^ol^i,'^, A^) 



The terms in Eq. (|43)l are ordered in such a manner that 
those in the left column can be later identified with terms 
bilinear in p"^ when using the GWT, whereas the terms 
in the right column are related to terms proportional to 
^10*^01 jg csgQi^tial for the following to note that 
terms with = /i or i/ = /I do not appear in the last line 
of Eq. (|43|l . i.e. they are zero. 

E. Energy kernel from the GWT 

The use of the SWT for the calculation of MR energy 
kernels relies on the construction of the very particular 
basis {di^,a^}. All practical applications, however, have 
used the GWT which typically provides, in any ar- 
bitrary single-particle basis, energy kernels of the form 
given by Eq. fTS]) . 

The next step is to compare the expressions of the 



energy kernel obtained from the SWT and the GWT. 
Using Eqs. and ([55]) . the transition density matrices 
are given by 

pOi = V yo* yo + V f/O- iii^iiil(44) 

V V \ U I / 

^ = E^t>^+E^^^-%§^(45) 

V y \ I / 

= E^^^+E^^^^-%§#(46) 

where the running index in the sums refers to the basis 
{a^,a+}. Inserting Eqs. (|44ll46p into Eq. Uni), the two- 
body part of the energy kernel is obtained from the GWT 
as 
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I'fJ, UfJL 

z^/^ \ I / \ I / 



9 Z^ >v¥'m0^'/'k ^I'P'^A'Ai /<I,„|*,\ /a,„|(f„\ ^zlZ^ ' 



(47) 



where each hnc of Eq. (|47p can be put in correspondence 
with those of Eq. (j43ll . Expression (|43|) is recovered term 
by term from Eq. ([T7|) as one can easily check that 

($o|$i,i/) ($o|$i,/i) ^ ($o|$i,;/, A^) 
($o|$i) ($o|$l) ($o|$i) ' ^ ^ 

for all pairs (/i, j/) such that v ^ [loxv ^ \1. However, one 
is left with additional terms in the last line of Eq. P7|) . 
i.e. for V = [i and v = pL, which have no correspondent 
in the last hue of Eq. (|4^ . 

In spite of such an apparent difference, the energy ker- 
nel Egwt[Q, 1] is strictly equal to Eswt[Q, 1] when calcu- 
lated from a genuine Hamilton operator, that is, within 
the strict projected-GCM approximation. To prove it, 
let us restart from last line of Eq. (|T7)) and isolate the 
possible source of discrepancies (i.e. terms with ly ~ fi 
and/or v = fl). Considering first the terms coming from 
^10*^01^ using the antisymmetry of the interaction and 
the properties of the matrix B, we find that 

iV- B* B* ($o|^i,i^) ($o|$i,^) 

^ Z^^v.v.'I'.'Pp ($o|$i) (4'o|«'i) 

+ 4Z^"v.v.0.0. ($o|«'i) (^ol^i) 

2^ ^.^.0.0. ($o|*i) ($o|$i) ' 

(49) 

which happens to cancel out exactly the contributions 
with [1 = 11 originating from p^^ p^^ . For this cance- 
lation to occur it is crucial that the same interaction 
is employed in both terms, as it is the case if a gen- 
uine Hamiltonian is used. Altogether only the term with 
V ~ p from the p*^^ p^^ contributions could not be com- 
bined with terms from k}^* k^^ . However, thanks to the 
antisymmetry of the two-body interaction, these terms 
also cancel out. 

Therefore, the prerequisites for a complete matching 
between estimates of the two-body energy kernel from the 
SWT and the GWT, independently of the appearance of 
divergences or not, arc that 



1. the same interaction kernel is used in both bilinear 
terms of the EDF in order to properly recombine 
specific terms coming from p*^^ p*^^ and k^^* 

2. the antisymmetry of the interaction kernel is prop- 
erly accounted for. 

The key point is that none of the two previous condi- 
tions are fulfilled in general when constructing a MR- 
EDF through an ansatz of the form of Eq. (fT5|) : the 
matrix elements w^^j,; are not necessarily antisymmetric, 
and v'^jj^i and v^Jl.^ are in general not the same. Conse- 
quently, the kernels £swt[0, 1] and £gwt[0, 1] defined in 
the EDF method by analogy with the Hamiltonian for- 
malism are different, with £gvi^t[0,1] containing terms 
that do not appear in £swt[0, !]■ 

V. CORRECTED ENERGY KERNEL 

The observations of the previous section are at the 
heart of the problems encountered in MR-EDF calcula- 
tions that rely on the GWT to construct the MR energy 
kernels from the SR-EDF. Let us analyze the situation 
further and demonstrate that there exist in fact several 
levels of problems of different physical origin. This will 
lead us to advocate the explicit removal of specific con- 
tributions to £gwt[Oi 1] in order to perform meaningful 
MR-EDF calculations. 

It is crucial to note already here that, as the correc- 
tion procedure proposed in the following is based on an 
analogy with the Hamiltonian formalism, it can only be 
implemented for functionals containing integer powers of 
the density matrices. We will come back to that crucial 
point when discussing the long term implications of the 
present work. 

A. Bilinear functional 

Let us start with the bilinear energy kernel £gwt[0, 1] 
given by Eq. (US)). As already discussed, generalizing the 
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functional to the multi-reference case using the GWT 
as a motivation leads to the same structure as the SR- 
EDF given in Eq. ([9]) , except that intrinsic densities have 
been replaced with transition ones. Again, the terms 
fg.^y[0,l] and SqIytI^A] might contain different ver- 
tices. 

We now introduce more compact expressions of the 
transition density matrices using the identity ($o|^i) = 



(50) 



where the sums run over all states. It is crucial to realize 
that the 8 terms in Eq. ()5m54p are arranged in exactly 
the same manner as the 8 terms in Eq. (|43|) . 



1. Correcting for divergences and finite steps 

By analogy with the strict projected-GCM case dis- 
cussed in Sec. Hvl the contribution to fGWT[0, 1] from 
the terms v = ^ and i/ = p, in the last line of Eq. ((M]) 
should be zero. However, the contribution encoded into 
the bilinear term p^^ p^^ is not zero for p, = v in line (j54p 
when the effective vertex v^p is not antisymmetric. Also, 
the vertices v^p and v'^'^ have no reason to be the same 
in an EDF approach; as a result, the contributions from 
pOipOi for ^ = jy and for {p = p, = D) do not 

cancel out anymore in line (|54p . 

This difficulty is caused by the use of the GWT as a 
motivation to construct MR energy kernel from the SR 
EDF. However, it must be made clear that the GWT 
itself cannot and should not be blamed for this failure, 
as it is stretched beyond its applicability when used to 
motivate the form of an energy kernel that does not cor- 
respond to the matrix element of an operator. It is prob- 
lematic though as the terms that do not cancel out in 
line ([M]) arc at the origins of the divergencies and steps 
seen in MR-EDF calculations based on the GWT ^ as 



where p and k are the intrinsic density matrices asso- 
ciated with |3>o}i whereas the matrix Z has non- vanishing 
matrix elements of the form Zp^, = {BjypA~^)* . Note that 
Z = DZD* = [BA-^Y is nothing but the Thouless ma- 
trix [23, H^l associated with the Bogoliubov transforma- 
tion connecting |$o) and Equation ([50]) is an alter- 
native form to the one given in Appendix E of Ref . (36l | . 



Using the mixed basis introduced in Eqs. (|40ll4ip . as 
well as Eqs. (|47l) and ([50]) . the interaction part of the 
energy kernel can be written as 



(51) 
(52) 
(53) 
(54) 

I 

discussed in detail in Ref. [2^ . 

As a matter of illustration, we now focus on diver- 
gences. Divergences will occur in £gvkt[0, 1] ($o|'I'i) if 
the two vacua involved are orthogonal; i.e. when it exists 
|$o) and |$i) in the MR set such that (^ol^"!) = 0. As 
can be seen from Eq. ((3T|) . the overlap (<I>o|$i) is zero if 
at least one of the matrix elements A^^ is zero. Because 
the non-zero contribution to fGH^Tp, 1] (<i>o|^i) coming 
from the term involving the pair (p, p) in Eq. (j54p is pro- 
portional to 

{Z,,f ($o|<I>i) = Coi {B;^f n ^I'p' ' (55) 

PP p'>0 

it will diverge as A*p goes to zero, except if another fac- 
tor A*,p, in the numerator happens to be zero as well. 
For particular MR calculations and very specific situa- 
tions, it is indeed possible that several conjugated pairs 
(P,P), {P',P'), ■■■are such that A*^ = A*,^, ==... = at 
the same time. In such a particular case, there will be 
no divergence of the MR-energy due to dangerous terms 
in Eq. ([5i)) . 

A regularized energy kernel is obtained by removing 
the spurious contributions to £q\yt ^^"-^ ^gwt that con- 
stitute the difference between £gwt[^, 1] and £swt[^, !]• 
This amounts to removing terms involving only one con- 



J 
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jugated pair (j),p) at a time in Eq. ((5^ 



PP , yPP 

(PpiPp4>p4>p ippipptpptpp 



p>0 

'^v^piPp<t>p4>p ~ ^v%p<t>pit>p) i^pp) ' (56) 



p>0 



(57) 



where the sum runs over p ^ i/ > only because all the 
matrix elements involving a given conjugated pair have 
been explicitly spelled out. 

Some important comments should be made 



(i) The sum of Eqs. (|56)) and ((57)) is of course zero if 
the MR energy kernel is obtained as the expectation 
value of the Hamiltonian, while it is non zero for a 
general functional energy kernel. 

(ii) The correction established through Eqs. ([5^ 
and ([57]) is zero when |$o) = Indeed, the Bo- 
goliubov transformation connecting the two states 
is trivial in this case with A being the unit ma- 
trix and B being zero. As a result, the correction 
does not modify the diagonal kernels £"[0, 0], which 
preserves the continuity between the SR functional 
and MR energy kernels S[p°°, = £[0, 0]. 



(iii) When pairing is not considered explicitly in the the- 
ory, the only terms to remove from Eq. (|54p are 
those obtained from p^^ p^^ for fx = u. As already 
discussed, such terms are indeed different from zero 
if the vertex v^p is not antisymmetric. As a matter 
of fact, we believe these terms to be responsible for 
the divergences seen recently in the restoration of 
angular momentum from cranked HF states [30l |: 
see Sec. IVIll 

(iv) The correction to £gwt[^^ 1] is independent of the 
normalization factor Coi and therefore of the phase 
conventions chosen to define the two vacua |$o) 
and |$i). This is so because the ratio Bpp/A*p is 
independent of Cqi. However, it should be kept in 
mind that the kernel £gwt[0, 1] is to be multiphed 
eventually by ($o|'J'i) where Cgi enters explicitly. 

(v) Equations. ([56]) and ((57|) not only remove possible 
poles leading to divergences but also correct the en- 
ergy kernel away from those potential poles. It is a 
crucial result of the present work to realize that cur- 
rent calculations are not only compromised through 
divergences but also through finite spurious contri- 
butions. 

(vi) As just said, poles associated to zeros of the norm 
overlaps are only a part of the pro blem. From that 
point of view, the study of Ref. [831 about the nodal 
lines of the overlap between cranked HFB states 
rotated by different values of the Euler angles is of 



prime interest. In Ref. |88| . it is shown in particular 
that the structure of nodal lines becomes richer as 
the states in the MR set are cranked to higher spins. 



2. Correcting for self-interaction processes 

In the last section, we have isolated the spurious con- 
tributions which are specific to £gwt[0, 1] when it is not 
obtained as the matrix element of the Hamiltonian. In 
particular, it was argued that the construction of MR en- 
ergy kernels based on the SWT is safe from divergences 
and steps. However, there exist others, less dangerous, 
spurious contributions which are common to £swt[0, 1] 
and £gwt[0, !]• This underlines the fact that the SWT 
is also stretched beyond its applicability when used to 
construct EDFs that do not correspond to the matrix 
element of a genuine operator. 

The first of those problems relates to the fact that 
the matrix elements v^^ might not be antisymmetrized. 
This leads to spurious self-interaction processes in the 
functional; i.e. the fact that a nucleon interacts with it- 
self [13 . The spurious terms not already removed by the 
correction discussed in the previous section, and which 
should be zero in £[0, 1] are obtained for v = fi ~ p or 
V ~ p = p in the terms that are independent of Z [line 
([ST])] or linear in Z [second and third line (|52ll53p ] 



9 / y \ fpVpVvVp fpVpVv'Pp 
p>0 



p>0 

IS 

p>0 



V-- -v'''' \ Z - 

9 / J \ fpipp't'p'Pp ifpippct>pipp j PP 
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9 z ^ \^ ^p^p^p4>p ^pipp^p(pp 



0j^PP-(58) 



It is clear that these terms are unphysical, as they vio- 
late Pauli's principle. This gives a strong motivation to 
remove them from the energy functional, as is sometimes 
done in DFT for electronic systems [1, Hi, H HH, llll , or 
to construct energy functionals that are self-interaction- 
free in the first place [11]. When using the standard 
correction, this will modify the SR energy density func- 
tional to what is called an " orbital-dependent functional" 
in electronic DFT [o^, [IHl , and will lead to equations of 
motion that are more difficult to solve numerically than 
the usual ones. From a phenomenological point of view, 
the merits of self-interaction corrected energy function- 
als for electronic DFT have not yet been fully clarified, 
as they improve some observables, but may degrade oth- 
ers when compared to uncorrected functionals (o^ . It 
has also been pointed out that correcting for the one- 
body self- interaction might not be sufficient, as there also 
might be also n-body self-interactions as well [ol, [o^ . 

In the nuclear context, the existence of spurious self- 
interactions in commonly used energy density functionals 
has been pointed out before [l|, [93, [111, but was never 
studied in any detail in the published literature. 
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3. Correcting for self-pairing processes 

The second category of problems common to 
Sswt[Qj 1] and £gwt[Q, 1] is less obvious to isolate. In 
fact, it differs from the two previous sources of problems 
in the sense that it is subject to interpretation. 

It is based on the observation that in the strict HFB or 
projectcd/GCM methods, where the same vertex enters 
^^''''[0, 1] and E'^'^[Q, 1], the terms corresponding to {ly = 
p; ^ = p) and {v = P'jU= p) recombine in Eqs. (|5m53p . 
As explained in Ref . [28| , such a recombination of terms 
can be interpreted as the fact that two nucleons in a pair 
of conjugated states cannot gain additional binding 
energy, as compared to the same EDF calculation with 
no explicit account of pairing, by scattering onto itself. 
Such a spurious self-pairing process should in principle be 



excluded from the functional, although the actual impact 
of these terms on observables is not clear yet. 

Compared to GWT-related and self-interaction prob- 
lems, removing the self-pairing energy does not amount 
to putting specific terms to zero. On the contrary, the 
reasoning is to allow certain terms to recombine and to 
provide a finite, non-zero contribution to the energy ker- 
nel. It is the remaining finite contribution which is sub- 
ject to interpretation. Here, we follow the argument from 
Ref. [15 mentioned above that takes the calculation with- 
out explicit treatment of pairing as a reference point. As 
a result, the advocated correction amounts to replacing 
specific matrix elements of w"'^ by the corresponding ones 
of vPf. Eventually, the correction for self-pairing amounts 
to the subtraction of 
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(59) 



r 



from the energy kernel, Eqs. (|51ll53p . The first two lines 
correspond to the spurious self-pairing energy on the 
single-reference level, while the last four lines represent 
the additional self-pairing energy contributing to multi- 
reference calculations. 



4- Chosen strategy 

The spurious terms causing divergences and steps, 
Eqs. (|56p and (|57p . are the consequence of a technical 
problem, i.e. using the GWT beyond its strict domain 
of validity. This could in principle be avoided by using 
directly the standard Wick theorem to define the MR en- 
ergy kernels. Eq. (j43ll . However, it is not clear at present 
if this would be numerically feasible. While the direct 
evaluation of Eq. might be too costly as it requires 
the explicit computation of the overlaps (4>o|$i, v, n) for 
all pairs of states {v,^), calculating the correction. ([TO]) 
and ([57]) . is not, and should become a standard proce- 
dure in the future. As we demonstrate in Ref. 12811 for 



pure PNR calculations, correcting SGWT[^^ 1] for diver- 
gences and finite steps indeed has a visible effect on the 
binding energy even when not accidentally hitting a di- 
vergence. The correction fluctuates on the order of 1 
MeV and leads to much smoother deformation energy 
curves. 

From a practitioner's point of view, it is indeed manda- 
tory to correct for divergences and finite steps as they 
seriously compromise, even inhibit, meaningful MR cal- 
culations. Correcting for self- interaction and self-pairing 
terms would also be desirable, either through explicit 
construction of self-interaction and self-pairing free en- 
ergy functionals, or through the correction of existing 
parameterizations of the functionals. Indeed, such spu- 
rious contributions to the energy functional violate the 
Pauli principle and hence arc always unphysical. 

However, correcting for self-interaction and self-pairing 
terms, Eqs. ((58)) and ((59l) . is a complex issue, both con- 
ceptually and computationally, i.e. it modifies the SR 
functional in such a way that the equations of motion 
are computationally more difficult to solve. However, 
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it is clear that the impact of sclf-intcraction and self- 
pairing on observables should be scrutinized. Existing 
studies in the context of HK-DFT for electronic systems 

Hi, S iH, S [9^ demonstrate that it is not an easy 

task and that it does not necessarily lead to an improve- 
ment of the performance of the functionals. As can be 
seen from Eqs. ([55]) and (|59p . sclf-intcraction and self- 
pairing arc a consequence of the violation of the Pauli 
principle within a pair of conjugated states. Even when 
removing both from the energy functional, there remain 
terms in the energy functional that break the particle ex- 
change (anti)symmetry between states that are not con- 
jugated. Eventually, it should not be forgotten that there 
is a price to pay for the enormous calculational simplifi- 
cation that EDF methods bring to the many-body prob- 
lem. Loosing antisymmetry on the level of the two-body 
density matrix is part of it, and how far it should be re- 
stored has to be determined as a compromise between 
the required precision and the necessary computational 
cost. 



5. Fixing the freedom in the regularization procedure 

In any regularization procedure, the quantity that re- 
sults from isolating and throwing away an infinity neces- 
sarily depends on explicit or implicit choices that deter- 
mine its finite value. In the present case, two elements 
are essential to fix such a freedom. 

First, and as explained above, we choose not to sub- 
tract finite self-interaction and self-pairing contributions 
which contaminate both SWT- and GWT-based expres- 
sions of MR EDF kernels. The second clement relates 
to the basis used to proceed to the regularization. Al- 
though the GWT can be applied in any quasi-particlc 
basis, which makes it so powerful, the application of the 
SWT necessitates to find a specific basis to express j^i) in 
terms of |<I>o), which we have managed to do as explained 
in section HVBI Of course, other bases exist that allow 
the application of the SWT. Onc^ such basis [2l,|7i| relics 
on diagonalizing the matrix Z^'^ rather than on apply- 
ing the BMZ decomposition to the Bogoliubov transfor- 
mation connecting |(f>o) and In the resulting basis, 
the comparison between SWT- and GWT-based formu- 
lae for operators matrix elements works just as explained 
in section IIV El in the sense that terms proportional to 
($o|'&i, i^)(*o|<I>i, i^)/(<l'o|*i>i)^ are present in the GWT- 
based formula but arc absent in the SWT-based one. Re- 
moving those terms from the GWT-based kernel would 
again eliminate divergences and steps. However, this 
would lead to a different regularized energy kernel as the 
factors and wave-functions weighting such terms are not 



the same as in the basis we prioritize^. Consequently, 
the choice of the basis used to proceed to the regular- 
ization impacts the final regularized kernel. Going from 
one basis to another amounts to reshuffling finite spuri- 
ous contributions between the different lines of Eqs. [43l 
orHZl 

Knowing that the regularization procedure is neces- 
sarily basis-dependent, just as standard sclf-intcraction 
correction methods in DFT are [8§| , the arguments that 
led us to prefer the basis we advocate are twofold. First, 
given |<I>o) and the application of the BMZ decom- 
position defines a unique basis, independently on the ac- 
tual representation of the two states. On the contrary, 
using the basis that diagonalizes Z^^ Z^ would require an 
extra argument to choose among the infinite number of 
vacua that can be used to express |$o) and Sec- 
ond, the basis we propose can always be found whereas 
it exists (rare) cases for which one cannot diagonalizc the 
matrix Z^^ Z^ [7|. 

Eventually, the motivated choice of the basis we ad- 
vocate and the decision to postpone to later the more 
challenging correction for self-interaction and self-pairing 
processes are the two elements used to fix the freedom 
that accompanies the present regularization procedure. 



B. Correcting higher-order functionals 

The cure proposed in the previous section to the prob- 
lems faced when constructing a MR energy density func- 
tional based on the GWT is specific to bilinear function- 
als. However, realistic functionals contain higher-order 
dependencies on the density matrices associated with 
many-body correlations and three-body forces. Those 
higher-order dependencies generate additional spurious 
contributions to the energy kernel £gwt[^: 1] which also 
have to be corrected for. The generalization of the strat- 
egy that we have followed is straightforward, as long as 
one considers integer powers of the density matrices. In- 
deed, one can formally relate them to higher-order multi- 
body forces in the Hamiltonian and identify, as we did 
for two-body forces, which terms are zero when using the 
SWT and replaced by non-zero terms when using the 
GWT. However, working out the SWT becomes lengthy 
and rather tedious for multi-body forces. 

Thankfully, one does not need to apply the SWT ex- 
plicitly, but rather proceed backwards using the con- 
nection between the SWT and the GWT discussed in 
Sec. IIV El This is more convenient because the GWT 
is easy to apply to any Hamiltonian. For a iV-body 
force, the SWT is recovered from the GWT by express- 
ing n-body transitions densities, with n S [0, A'^], through 



* In fact, there is an infinity of such bases associated with the 
freedom to choose the third vacua with respect to which | #0) and 
l^i) must be expressed when using the technics of Refs. [23I. |79|| . 



One noticeable exception is PNR for which all convenient bases 
are the same. 
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Eqs. (|44II46P and by using the identity 

($o|$i) (*o|<i>i) ($o|$i> ^ 

which gcnerahzes in passing the notations introduced in 
Eqs. p2j|33p in an obvious manner and which is vaUd as 
long as all f^'s belong to different conjugated pairs. When 
two or more i^^'s belong to the same conjugated pair, the 
validity of Eq. [SD] is lost and the left-hand side factor 
should be replaced by ($o|'l'i, J^i, • • • i^n) = 0. Using 
such shortcuts, one can obtain an explicit "SWT" form 
for the matrix clement of multi-body forces in a very 
economical way. 

Let us illustrate the above strategy for a three- 
body interaction. This will correspond to complement- 
ing the SR-EDF and the energy kernel 
Egwt[^, 1] with a trilinear component. The three-body 
interaction is written as 

wi23 = X! ^y'^'™™ ajaljolanCiraai , (61) 

ijklmn 

where the three-body force matrix elements are fully an- 
tisymmetric. Within the strict projected-GCM approxi- 



mation, the application of the GWT leads to 

— L \^ T, riOl 

ijklmn 

+ i ^ijklmn P^i tjfc «^mri ■ (^2) 

ijklmn 

Operating as described above leads to the expression of 
the overlap (see Appendix [X]) as it would be obtained 
from Eq. ([^^ using the SWT. Such an expression com- 
plements the last line of Eq. for the two-body Hamil- 
tonian. 

As before, we now define the trilinear energy kernel 
^gw'tI'^' 1] and Sq^j^tI^' 1] within the EDF approach in 
close analogy to the two terms of Eq. (|62p . The two con- 
tributions are expressed in terms of the effective vertices 
yppp and ■yP'*'', respectively. Just as for the bilinear parts, 
^givtI^: 1] and £^.'^^[0, 1] need now to be corrected for 
spurious contributions. Proceeding in a similar way as for 
the bilinear functional, and using the SWT as a reference 
point, Eq. (jAl|l . the terms to be removed are 
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[ZppY Z-,,. (64) 



We have not made use of any antisymmetry properties of the vertices in Eqs. (|63p and (j64p . even though the vertex 
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^PKK ^g^^ always be chosen to be antisymmetric, at least 
with respect to the second and third indices on the one 
hand and to the fifth and sixth indices on the other. 

A few observations similar to those made for the cor- 
rection of the bilinear functional can be made for the 
trilinear functional. First, ^^^[0, 1] and £^^q[Q,1] sum 



up to zero if, and only if, both effective vertices refer 
to the same fully antisymmetric three-body interaction. 
Second, both correction terms are zero for [$o} = I'i'i) 
and do not modify the underlying diagonal kernel £[0, 0], 
keeping valid the link between the SR functional and the 
MR one £[p, k, k*] = £[0, 0]. Third, the correction to the 
functional kernel £gvkt[0, 1] is independent of the nor- 
malization coefficient Coi- 

Note finally that the three-body terms should also 
be corrected for spurious self-interaction and self-pairing 
problems. Expressions equivalent to ([58]) and ((59l) can 
be derived without difficulty for three-body vertices (not 
shown here). 



VI. APPLICATION TO PARTICLE NUMBER 
RESTORATION 

In the present section, we specify the previous findings 
to PNR calculations. A more extensive discussion of that 
particular case, including results of realistic calculations, 
is presented in Ref. [1^. When considering one kind of 
particles only, the MR set appropriate to PNR is given by 
the ensemble of quasi-particle vacua {\^^);(p £ [0, 27r]} 
which correspond to states rotated in gauge space by an 
angle ip. The MR energy functional that amounts to 
restoring the particle number N is given by 



it is worth pointing out that Eq. ([67|) corresponds to a 
state |$o) = Co Hi/ Q^i^lO); where the quasi-particle oper- 
ators {q!i/, aj} arc defined through a BCS-type transfor- 
mation 



with V = p or p and where we have defined 
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The second vacuum required to construct the energy 
kernel £[0,(p], i.e. where here |$i) = 1$,^). This state is 
obtained by rotating |<i>o) in gauge space by an angle if, 
i.e. 



1$. 



5^'^*|$o) = Yliup + vpc^'^'^a+a+m. (69) 

p>0 



From a technical point of view, particle-number restora- 
tion is simple because all vacua belonging to the MR set 
share the same canonical basis. As a result, the Bogoli- 
ubov transformation linking any pair of vacua in the set 
is itself canonical, i.e. 
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It is also interesting to compute the transition densities 
through Eq. ([50)1 . For example, it leads for the diagonal 
normal transition density matrix 
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which is real. 

In the previous equations, the state |$o) is the quasi- 
particle vacuum at gauge angle ip = 0. Focusing on 
projection after variation, such a state is obtained from 
a self-consistent (possibly constrained) SR calculation. 
Considering the case of an even-even nucleus, the vac- 
uum state, possibly breaking time-reversal invariance, is 
written in its canonical basis {a^, af} as 

\'^o) = Y[{up + Vpa+a+)\0), (67) 

p>0 

where [0) is the particle vacuum and where {up;Vp} are 
real numbers. According to our convention, the product 
in Eq. ([67]) only runs over the "positive" half of the basis. 
The state |$o} is normalized, with the convention that 
Up + Vp = 1. To underline the link with previous sections. 



Ov ^ yO_* yO_ ^JjO Em V"- 

I pp pp pp ' pp ^ ^ pp 



"p 



vl e^^v ' 



(73) 



whereas the anomalous ones are obtained in the same 
way as 
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Before turning to the correction, let us discuss the rel- 
ative phase between the two vacua. On the one hand, 
applying Eq. ([31]) provides 
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whereas on the other hand, the explicit canonical 
forms given by Eqs. (|67ll69p lead to ((f>o|(f>ip) = 
Ylp^o {ul + vle"^^^) . Thus, the phase convention asso- 
ciated with the canonical forms is equivalent to having 
chosen Cqi Hp^o e"**^ = 1. 

The fact that the Bogoliubov transformation {A, B) 
is canonical from the outset is a specificity of PNR. In 
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particular, it corresponds to having the trivial unitary 
transformations C = D = 1 in Eq. (QH). For almost 
all other cases of interest, however, i.e. Angular Momen- 
tum Projection and/or GCM-type Multi- Reference cal- 
culations, A and B cannot be obtained analytically. This 
means that the matrices A and B must be computed us- 
ing Eq. for each pair of vacua belonging to the MR 
set. To this end, methods that will be outlined in Scc. lVIII 
have to be applied. Their implementation is underway 
and will be discussed elsewhere. 

Finally, the upper and lower components associated 
with the quasi-particle operators {ap,a^} introduced 
through Eqs. (|40II41|) take the form 



~Vp \p) , 
Up \p) , 
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Vp \p) , 
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where \p) and \p) denote the single-particle states created 
respectively by a+ and a^, respectively. 



1. Bilinear functional 



The different types of spurious contributions to the 
bilinear energy kernel read in the PNR case as 
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Subtracting £^q[0,(^] and £qq[Q,ip] from the PNR energy kernel does not call for a modification of the imderlying 
SR functional whereas removing 1^5/ [0, ip] and fgpp, ip] would. 



2. Trilinear functional 

The spurious contribution to be removed from the GWT-based EDF kernel £ppp[0, ip, Lp'\ is given by 
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In this expression appears a transition density matrix as 
defined by Eq. ([73]) and evaluated at a second gauge an- 
gle ip' . This second gauge angle refers to the fact that the 
isospin degree of freedom must be treated explicitly when 
correcting the trilinear functional, as there will be terms 
that are bilinear in one isospin and linear in the other. As 
a result, both projections on neutron and proton num- 
bers must be considered explicitly. The isospin could 
be omitted for the bilinear functional because, as long 
as we do not mix neutron and proton in the mean-field 



and only deal with neutron-neutron and proton-proton 
pairing, the conjugated states p and p always refer to 
the same isospin and are rotated by the same angle Lp 
in gauge space. In Eppp and f'"', however, the particle 
V may have a different isospin from the one of the pair 
{p,p) and a second gauge angle ip' must be attached to 
it. 

Finally, the contribution to be removed from the func- 
tional energy kernel f '"^ reads as 
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where, again, ip' ^ li and only if v and p have the same 
isospin; otherwise Lp' ^ Lp. 

VII. NO-PAIRING CASE 

Multi-reference EDF calculations are sometimes per- 
formed without an explicit account for pairing correla- 
tions; i.e. the reference states of the MR set take the form 
of Slater determinants rather than quasi-particle vacua. 
We concentrate on such a case in the present section. 
This is of particular interest since recent MR calcula- 
tions based on triaxial cranked Slater determinants and 
aiming at restoring angular momentum displayed diver- 
gences [1^ which we believe to be related to the problems 
discussed in the present work. 

To investigate the zero-pairing realization of the for- 
malism outlined above we consider a MR energy kernel 
associated with two Slater determinants |$o) and |$i} 
given by 

AT 

i*o) = n«^io) ' 

1=1 

N 
i=l 

where N represents the number of particles whereas 
{a^}i=i.Ar and are the creation operators as- 

sociated with occupied, i.e. hole, single-particle states 
{\ai)}i=i,N and {|&i}}i=i,Ar, respectively. In the general 
case, we have expressed the two quasi-particle vacua, 
Eq. (PP]) . in a common, arbitrary single-particle basis. 
However, numerical applications often make use of two 
non-equivalent sets of single-particle states to define the 
two vacua. This is naturally the case when dealing with 
Slater determinants. Thus, and as a by-product, the fol- 
lowing discussion of the zero-pairing case will sketch how 
to proceed when the quasi-particle vacua are defined with 
respect to different single-particle bases. 

An important result of the present work (Section IIV[) 
consists in finding a convenient basis through the BMZ 
decomposition of the Bogoliubov transformation con- 
necting two quasi-particle vacua, Eq. (P7|) . The Slater 
determinant limit discussed in the present section offers 
a chance to illustrate more intuitively how this works. 



The single-particle subspace spanned by the 
{|&i)}i=i,7V is a priori different from the one spanned by 
the {\ai)}i=i^N ■ However, it is clear that the subspace 
spanned by the addition of the two bases is at maximum 
of dimension 2N . Therefore, one can always complete 
the set of hole states of |$o) by an appropriate choice 
of N particle states, denoted by {\o.^)}k=i,N , in such a 
way that 

N N 

|6,)==^|a,)(a,|M+^|as)(afe|fe,). (85) 
i=i fc=i 

Similarly, one can introduce a set of TV particle states for 
|$i), denoted by {\bk)}k=i,N such that 

N N 

\a^)^Y.\^J){h^a,)+Y.\h){h\<^^)■ (86) 
j=i fe=i 

In practice, the most convenient choice is to associate to 
each state \bi) an intermediate state \di) given by 

N 

\d,)^\h)-Y.\a,){a,\h). (87) 

Then, the N states \aj.) are obtained by a Schmidt or- 
thogonalization of the N states \di). If the state \hi) can 
be completely described by the hole states of |<&o)j the 
associated \di) cancels out and the dimensionality of the 
considered matrix can be reduced. 

We associate the creation operators and 6^ to the 
states I a J,) and respectively. One has in particular 

|$o) = and |$i) = 0, for aU A: = 1, . . . , iV. 

Starting from the truncated space as introduced above, 
a formalism adapted to MR- EDF calculations with Slater 
determinants could be introduced without referring to 
quasi particles at all. Then, expressions of the functional 
kernels and the corrections to spurious processes could 
be obtained directly in the language of particles. How- 
ever, our present goal is to demonstrate how the formal- 
ism written previously in terms of quasi particles can be 
directly applied to the particular case where pairing is 
not considered explicitly. The first step is thus to match 
the particle creation/annihilation operators with quasi- 
particle creation/annihilation operators for Slater deter- 
minants. Although this can be found in textbooks [36| . 
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we present it below for the sake of a self-contained de- 
scription of our method. 



A. From particles to quasi particles 

The matching appropriate to the notation used in 
previous sections is achieved using the standard quasi- 
particle representation of Slater determinants. We in- 
troduce the two sets of quasi-particle annihilation oper- 
ators ai, and /S^ and restart from the general form of the 
transformations between particles and quasi-particles op- 
erators: 




4N 




(88) 



and 




(89) 



Here, the index AN is used to stress that only a specific 
finite number of quasi-particle states is necessary. The 
annihilation operators verify a^, |$o) = and |$i) = 
for 1/ = 1, . . . , 2N. It is clear that this relation is fulfilled 
if ai, (resp. (3i,) is proportional to a linear combination 
of hole creation operators of |$o) (resp. I^*!)) or to a 
combination of its particle annihilation operators. This 
leads to a simple block structure^ for the matrices C/°/^ 
and 




■2N 




(90) 



2N 



which is obtained thanks to a specific ordering of hole 
and particle states, that is, we have {a'^}i/=i,2N = 
({a+}i=i,Ar,{a^}^=i,Ar)- 



Bogoliubov transformation between two Slater 
determinants 



To express the matrices A and B associated with the 
transformation between the two sets of quasi-particles 
{a,y,al} and {/3^,/?^} introduced above, we first define 
the transformation 





(91) 



It is worth mentioning that there is a flexibility in choosing the 
lower-right block and upper-left blocks of [/"/l and V"/! respec- 
tively. Here we choose the simplest prescription which in our 
opinion is also the most convenient in applications. 



where R is the overlap matrix that can be schematically 
written as 



R 



m\a,)} {{h\aj)} 
{{b,\a,)} {{b,\aj)} 



2N 




The matrices A and B can finally be deduced using 
Eqs. (Ell), (IMl) and (Ell) 



(93) 



A = U^'+R+U^ + V^+R^V\ 
B = V^^ R+U^ + U^^R^V^ , 



which, using the simple form of U'^^^ and leads to 

■ (94) 



A 




B 



2N 




2N 



C. Bloch-Messiah-Zumino decomposition 

The matrices A and B fulfill all usual relations asso- 
ciated with Bogoliubov transformations. Guided by the 
standard BMZ theorem, we introduce the two matrices 
S and T defined as 



S 
T 



B*B'^ . 
-AB+ 



B*A' 



(95) 



These two matrices play a role similar to that of the 
normal and anomalous density matrices in the SR case. 
In particular, we have the following relationships 



S^-S = -TT^ 
ST = TS* . 



(96) 



However, it should be noted that the matrices S and T 
are not related to the transition density matrices in any 
trivial way. Using the expression of A and B, S can be 
written as 



while T takes the form 





-{AB+f 



2N 




,(97) 



2N 



AB+ 












/ 2N ^ 




:) 



2JV 



in such a way that T is antisymmetric as expected. In 
the following, it will be useful to have a more explicit 
form of the matrices , and T. Using Eqs. ((92l) 
and (|85ll86p we deduce 



S, 



+ (a.|(l-/^)|a,r , 



4 = 



(a.|p"|a,)' 



(99) 
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One can see from Eq. ([99)) that all the information rel- 
ative to the overlaps between the single-particle wave- 
functions of the two bases is encoded into the three ma- 
trices S^, and T. 

The first step to obtain a BMZ decomposition of A and 
B is to diagonalizc the matrix S. The block diagonal form 
of S, Eq. ([57]) . shows that the unitary transformation 
allowing to do so is also block diagonal; i.e. diagonalizing 
S is achieved by diagonalizing in the subspace of hole 
states of |$o) and in the subspace of particle states of 

$o)- As a result of such a transformation, new hole and 
particle states are obtained for |$o)i denoted by \ai) and 

flj), respectively. We define the corresponding creation 
and annihilation operators by df and a, as well as the 
eigenvalues of S-^ and by Af and A^, respectively. 
Expressing the two relationships of Eq. ([M]) in the basis 
where S is diagonal leads to 



N 



4(4-1) 



Ar(Ar - 1) 



^j(Af-A^) 



i=l 
N 



2 



2 



(100) 

(101) 
(102) 



which shows that %j may differ from zero only if Af = Aj . 
Let us now consider the different possible eigenvalues of 
. There exist three different cases 

• Af = 0: this means that the cigcnstatc \ai) of S-^ is 
orthogonal to all particle states of |$i) and there- 
fore can be written as a linear combination of the 
occupied states in |$i). Eq. (|10ip automatically 
implies that T^j = for all 7 = 1, . . . , iV. 

• Af = 1: this means that the eigenstate \di) is fully 
contained in the subspace of particle states of |$i). 
Again, Eq. (jlOip implies that %j — for all j — 
1,...,7V. 

• 1 < Af < 0: the corresponding eigenstate is neither 
entirely contained in the space of particle states 
of l^i) nor in the space of its hole states. This 
also means that at least one matrix element Tlj is 
non vanishing, for j = I, . . . , N. Accordingly, there 
exists at least one eigenvalue Aj such that Aj = Af . 

The same classification could be made for the eigenvalues 
Aj except that Af = or Af = 1 correspond in this case 
to eigenstates |a,) belonging to particle or hole states of 
$1), respectively. 

It could be checked that the diagonalization of S is 
not affected by a unitary transformation that acts sep- 
arately among the particle or hole states of |$i). In- 
deed, p^^ and (1 — p^^) entering Eq. ([99]) are invari- 
ant under such a unitary transformation, respectively. 
To advance towards the BMZ decomposition the Bogoli- 
ubov transformation (|93|) . one now needs to transform 



the (quasi-)particle states associated with |$i). This can 
be achieved by repeating the above procedure for the new 
matrices S" = B^B* and T' = B^A 



5' 

and 
T' 



X^X* 



2N 








S' 



X 



, (103) 



2JV 



B^y 

-{B'^yV 



2JV 



T' 



I (1Q4) 

2N 



where the matrix elements now read 



c'^ 

T'- = 



+ (6,|(l-p°")|fe,) 

+{b,\p''\bjy , 



(105) 



The unitary transformation that diagonalizes S' consti- 
tutes the second step of the procedure. As for S, the 
transformation is block diagonal, in such a way that hole 
and particle states transform among themselves. This 
leads to a new set of single-particle (eigcn)states denoted 
by \bi) and |6j). The corresponding eigenvalues of S" 
are given by A'f and A'". A gain, we can classify single- 
particle states in three categories according to the values 
of X'^ for hole states and to the values of A'^ for particle 
states. 

From a practical point of view, no numerical diago- 
nalization of S' is actually needed. Starting from the 
properties of the states \di) and \dj) obtained through 
the diagonalization of S 



(a.|(l-p")|a,) = Af5. 



(106) 



one can deduce the hole and particle states of |$i) that 
diagonalize S'^ and S'^ . We consider here the case 
where the eigenvalues Af and Af differ from zero^. With 
each hole state \di), resp. particle state |a^), we asso- 
ciate a new particle state |&j), resp. hole state \bi), of 
$1) through 



1 ^ 



1-p 



11 



fc=i 

N 



\d,), (107) 



1%) . (108) 



Those states are 



If the eigenvalue Af/A^ is or 1, the state |ai)/|ai) is already 
a particle or a hole state of I"!?!). These cases can be treated 
through a preliminary step where some of the states jfe;) and \b^) 
arc directly identified as particle and hole states of |#o)- The 
other particle and hole states of are changed accordingly to 
fulfill orthogonality conditions. 
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• orthonormal and have overlaps with the states di- 
agonahzing 5* given by 



(AD 



1/2 



eigenstates of S'^ and S'^ and fulfill 

f! — f- 



(109) 



(110) 



From the above relation, we deduce in particular 
that Af = and Af = AJ. 

Now, a hole state with eigenvalues Af can now be 
decomposed as 

= |a,)(a,|6,)+ ^ (HI) 

where only the particle state joj) from which \hi) is con- 
structed appears, while the summation over holes of |3>o} 
is restricted to states that have the same eigenvalue as 
\hi). As a result, the number of states in the expansion 
is greatly reduced as compared to Eq. ([55]) . 

The matrix A is not yet written in a canonical form. 
To complete the BMZ decomposition, we perform a third 
transformation consisting of mixing the states \ai) (resp. 
fij}) within each degenerate subspacc of (resp. 5^) 
in such a way that the matrix A (resp. y) becomes diag- 
onal. The product of the transformations diagonalizing 5* 
and A provides the unitary transformation D introduced 
in Sec. IIVI where we explained our method in the general 
case. As A is being diagonalized, the states (resp. 
|6i)) are modified accordingly through Eq. (|107p (resp. 
Eq. (jlOSp )^. Together with the transformation that di- 
agonalized S", this provides the unitary transformation 
C introduced in Sec. HVl 

Finally, the BMZ decomposition of the Bogoliubov 
transformation between two Slater determinants (j93[) is 
achieved in the sense that the matrices 




B 



2N 




(112) 



2N 



Before 
Bloch-Messiah-Zumino 



After 

Bloch-Messiah-Zumino 



: 



4 



bi 



l*o) 




FIG. 1: (Color online) Schematic illustration of the different 
bases introduced in the text. Left: the Slater determinants 
|<I>o) and |"I>i) are expressed in their respective natural bases 
{|ai)} and {|&;)}. The shaded areas indicate that one state 
of the basis {!&;)} ({I'^i)}) a priori spread over all particle 
and hole states of the basis {|ai)} ({!&;)})• Right: the BMZ 
decomposition leads to two bases which are such that each 
transformed state \hi) (|5i)) decomposes only onto two states 
{ai,at) {{bi,bi)) (one hole and one particle). At the same 
time, the latter decomposition highlights that the notion of 
conjugated pairs for Slater determinants relates to the associ- 
ation of each hole state with a particular particle state. Note 
that energy levels have no specific meaning in this figure; they 
are just meant to characterize occupied and empty levels for 
each of the Slater determinants. 



are in canonical forms when expressed in the new bases. 
In the present case, conjugated pairs are made of a parti- 
cle state and a hole state af of |$o). Through a BCS- 
likc transformation, the latter pair is associated with a 
unique particlc-holc pair (b^b^) of 



b+ = Xs~a++%a+, 



(113) 
(114) 



where the matrix elements of A and B are given by the 
canonical expression 



Aji = 5jilJ)i\a.i) , Bji = 5Yi(bi\ai) 
= 5yi{a.i\bri) , yji = 5ji{ari\bi) 



(115) 



Finally, Fig. [T] illustrates the procedure and the meaning 
of putting the Bogoliubov transformation of dimension 
2A^ that connects two Slater determinants into a canon- 
ical form. 



D. Correction for spurious processes 

Given the BMZ decomposition of the transformation 
between two Slater determinants obtained above, we can 



Although the last step bringing A into a canonical form in- 
volves an additional transformation, wc use the same notation 
{|ai); \ai)}i—i^...^N and |t'!:}}i=i,...,jv to denote the bases 

before and after that last transformation. 



The bases {ai,a|} and arc the analog of the bases 

{oijj, qJ,} and {/3fi,/3^} defined by Eq. I I27I I in the general case. 
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apply directly the results derived in Sec. |TVl using the 
new basis {5,^, a|}, to write the state |$i) under the form 

N 

\<i>i)^l[{A*,+8:,~arat)\%), (116) 

and the transition density matrix between the two Slater 
determinants as 

N 

p"'=p"° + Y.\&,)Zu{a,\, (117) 
1=1 

where, as before, = B*^ A*~^ . As can be seen from 
Eq. (jll6p . the picture emerging in the new basis is very 
intuitive since the Slater determinant is obtained as 
a linear combination of up to. very specific, N particle- 
N hole configurations on the vacuum |$o). In partic- 
ular, each hole di is combined with a specific particle 
di] the latter association defining the notion of conju- 
gated pair for the Bogoliubov transformation that links 
two Slater determinants. Such an expression simplifies 
tremendously the general expression connecting two non- 
orthogonal Slater determinants where each hole state can 
be excited into any particle states, up to infinite energy. 

As for calculations with an explicit treatment of pair- 
ing, Eq. pi7p can be used to express the MR energy 
kernel on the basis of the GWT. This leads to 

N 

SgWt[0, 1] = ^ (taidi + td-ai 




+ 2 E ^Za,,.a,ZuZj,. (118) 

Starting from this expression, the self-interaction 
f^j[0, 1] and the spurious contribution to the MR en- 
ergy due to the construction of the energy kernels on the 
basis of the GWT £^%[0, 1] can be identified 

1 ^ 

1 ^ 
1 ^ 

S'c'gIOA] = 2E^ar<^.a,a,^S, (120) 

whereas, of course, no spurious self-pairing occurs in the 
present application. 



Again, both terms are zero for energy kernels obtained 
as the matrix element of a Hamiltonian. The latter of 
the two terms is very likely to be responsible for the dif- 
ficulties recently encountered in calculations aiming at 
restoring angular momentum using a MR set of cranked 
Slater determinants [s^]. We thus advocate the removal 
of the contribution given by Eq. (|120[) in such a context. 



VIII. CONCLUSIONS 

The present work concentrates on multi-reference 
(MR) calculations, customarily called "beyond-mean- 
field" calculations in the literature, performed within the 
Energy Density Functional (EDF) formalism. The multi- 
reference method is nowadays routinely applied with the 
aim of including long-range correlations associated with 
large-amplitude collective motions which are difficult to 
incorporate in a more traditional single-reference (SR), 
i.e. "mean-field", EDF formalism [l|. So far, the frame- 
work for such MR-EDF calculations was set-up by anal- 
ogy with projection techniques and the Generator Coor- 
dinate Method (GCM), that has been rigorously formu- 
lated so far only within a Hamiltonian/ wave-function- 
based framework [ssi . [36| . 

The first achievement of the present work is to demon- 
strate that the usual extension of the single-reference en- 
ergy functional k, k*] into the non-diagonal energy 
kernel £^[0, 1] at play in MR calculations through £[0, 1] = 
£[p°^,K°^,K^°*] is ill defined. The latter extension, 
based on the Generalized Wick Theorem (GWT) is 
well defined within a Hamiltonian/wave- function-based 
projected-GCM framework, but happens to be at the 
origin of spurious divergences and steps in MR-EDF cal- 
culations, as recently realized for particle number and 
angular momentum restorations [2^, [2^, [s^l • 

The second achievement of the present paper is to pro- 
pose a method to identify, for any type of symmetry 
restoration and/or GCM-based calculations, the spuri- 
ous terms in the MR-EDF responsible for divergences, 
steps, self-interaction and self-pairing. The versatility of 
the method also allows to take care of difficulties encoun- 
tered when mixing within the MR-EDF calculation states 
obtained from quasi-particle excitations [3l[. 

In practice, the method requires first to put the Bo- 
goliubov transformation connecting two vacua |3>o) and 
|3>i) into a canonical form. In the corresponding canon- 
ical basis, spurious contributions to the MR energy ker- 
nel £gwt[0,M can be identified and, if necessary, re- 
moved. Subtracting the terms responsible for divergences 
and steps is a prerequisite to perform well-defined multi- 
reference EDF calculations, and that is what we advocate 
in the present work. However, the quantitative impact of 
self-interaction and self-pairing processes on observable 
as well as the technical difficulties associated with their 
removal remain to be studied. As they do not lead to 
dramatic consequences, i.e. divergences and steps, this 
can be postponed. 
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As the proposed procedure must be implemented for 
all pairs of vacua involved in the multi-reference calcula- 
tion, the correction increases the basic computation cost. 
On the other hand, the number of summations in the cal- 
culation of the MR energy is significantly reduced in the 
canonical basis extracted to remove the spurious contri- 
butions. Thus, it might be of interest to take advantage 
of that fact in actual implementations of the MR-EDF 
formalism. This is the goal of a forthcoming publication 
to discuss such an implementation. In the meantime, an- 
other publication is devoted to illustrating the features 
of the correction method proposed in the present paper 
by applying it to the rather transparent case of Particle 
Number Restoration [H]. 

Despite the success of reg ularized MR-EDF calcula- 
tions documented in Ref. [28| . some important questions 
remain unanswered as of today regarding the foundation 
and the implementation of MR methods within a con- 
sistent EDF framework. First, it is unclear in which 
sense the "projected" functional that one manipulates 
in such calculations can be attributed, at least implicitly, 
to a state belonging to a specific irreducible represen- 
tation of the symmetry group of the nuclear Hamilto- 
nian. The latter statement may seem rather paradoxical 
as the main purpose of those calculations is precisely to 
"restore the symmetries". The crucial point is that, as 
opposed to the strict Projected Hartree-Fock-Bogoliubov 
approach [33,[3^[3^, the MR energy functional is not ob- 
tained as the average value of a Hamiltonian in the pro- 
jected state but as a functional of (non-observable) tran- 
sition densities defined for each pair of vacua belonging 
to the MR set. The rather complicated dependence on 
those densities and the use of different effective vertices 
(interactions) in different channels of the functional for- 
bid the explicit re-factorization of the energy in terms of 
the projected state. The benefit of such a method is that 
correlations which go beyond the strict projected-GCM 
approximation can be easily incorporated into the MR 
formalism. Again, the drawback is that it is not clear in 
which sense symmetries arc actually restored in existing 
MR-EDF calculations. The clarification of such a ques- 
tion and the proper formulation of symmetry restoration 
within a well-defined MR-EDF theory is mandatory in 
the near future. 

Last but not least, the present work raises three impor- 
tant questions; the first two being related to the construc- 
tion of energy density functionals; i.e. Skyrme, Gogny or 
any other types of realistic functionals. 

(i) Only the removal of the most dangerous spurious 
terms associated with the use of the GWT as a ba- 
sis to construct phcnomcnological MR energy ker- 
nels has been advocated in the present paper and 
implemented in Ref. [28| . However, it is worth 
noticing that SR- and MR-EDF calculations are 
also plagued with less dangerous pathologies re- 
lated to self -interaction [sj and self-pairing [2^ . 
In particular, such pathologies contaminate MR en- 
ergy kernels independently on whether the SWT or 



the GWT is used as a basis to define them. At 
this point though, we do not advocate the correct 
for self-interaction and self-pairing processes briefiy 
discussed in the present paper. The reason is that 
such additional corrections will impact explicitly 
the SR functional on which the MR calculation is 
based. This means that self-interaction- and sclf- 
pairing-free SR functionals must be constructed in 
order to test the importance of such spurious pro- 
cesses. This is a non-trivial step to take because 
doing so will modify the structure of usual function- 
als and most importantly the way self-consistent 
SR calculations usually work [s^. Such a work is 
underway. 

(ii) The correction method proposed in the present 
work is based on using the Hamiltonian/ wave- 
function framework, together with the Standard 
Wick Theorem [4l|, as a reference point. A non 
trivial implication is that the procedure provides a 
way to correct functionals which depends only on 
integer powers of the normal and anomalous den- 
sity matrices [1^. As we do not see at this point 
how to proceed otherwise, the present work acts 
as a strong motivation to construct SR function- 
als that only involve integer powers of the density 
matrices in the near future. For reasons that are 
well known to practitioners, functionals depending 
on the third power in the local density do not work 
well enough and it remains to be seen whether or 
not using fourth powers of the local density is suf- 
ficient to construct high-precision energy density 
functionals. 

(iii) The present work is a satisfactory solution to an 
extreme problem faced by MR energy density func- 
tional calculations. However, it is not entirely 
satisfying from a fundamental point of view since 
it amounts to correcting a phcnomcnological con- 
struction of multi-reference energy kernels that is 
ill-defined in the first place. Combined with re- 
lated issues regarding restoration of symmetries, 
the present work calls for formal developments to- 
wards the derivation of a MR-EDF formalism from 
first principles. In particular, the regularization of 
the energy kernels happens to be basis-dependent, 
just as standard self-interaction correction meth- 
ods in DFT are [s^. A nice feature of the practi- 
cal method we use to proceed to the regularization 
is that it leads to a unique basis among all pos- 
sible ones. It is so because, given |$o) and l^i), 
the application of the BMZ decomposition defines a 
unique basis, independently on the actual represen- 
tation of the two states. Still, the fact that the re- 
moval of spurious contributions is basis-dependent 
in the first place is not satisfactory and deserves 
some further attention in the future. 
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APPENDIX A: EXPRESSION OF ($011)1231*1) 

Using the method depicted in Sec. IV B| we obtained 
the following expression: 



+J E u^Mv^^aiylVuxV^x s;,($ol<f>i, 

+\ E yil*y^u°rn-^y!^.ui^ykx B^,s*,,(ci>oi$i, a) 

+\ E u^MK-;^!^ulxykx ^fcimn 5*,%B^,($„i$i, z., A) 
+\ E yll*y^ly'M*A^nx ($oI$i) 
+i E yil*y'.y!Ml^ixUnx B|,(<i>oi<i>i, a) 
+i E yll*y^ly!^yk-^yn:xuL s*^($oI$i, 
+i E u^Mv°M;v^\uL %fei™n s*,($oi<f>i, z.) 

+ i E yil^y^y^^S'k-^U'^nxUnX V^.Umn ($0 |$i , ^i. A) 

+ i E aC^°A %feimn S;,B|^($o|<i>l, i^. A) 

+ i E ^M^fcpK^Af>°A %M™„ B;,%(<i>o|<i>l, fi) 

+ \ E ^.°^M^°m^™aC^°A ^'^Jklmn S*,%B|,(<i>o|$l, i^, A^, A) (Al) 



where the sums run over all the indices, with the con- 
vention that when two or more t'i's belong to the same 
conjugated pair, ($o|$ii ^^ij • • • ~ 0. In the latter 
expression, the first 8 terms come from the expansion 
of p^^ p^^ p^^ . Among these terms, the four last one will 



provide spurious contributions to £q]^t- Similarly, the 
remaining 8 terms originate from p^^k}'^ k'^^ whereas the 
last 4 of them provide spurious contribution to S^^j,. 
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